Objectifs de la formation

Note

La présentation est disponible à l’adresse: https://taq.info/trainingMetabarcoding/

  1. Comprendre le séquençage Illumina et les fichiers FASTQ
  2. Manipuler des données de séquençage avec R
  3. Découvrir BARQUE, un pipeline de métabarcoding québécois
  4. Appliquer les concepts sur des exemples pratiques

Avant de commencer

Installation des librairies R requises

  1. Ouvrir RStudio
  2. Envoyer les commandes suivantes dans la console
# Installation des packages Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("ShortRead", "Biostrings", "Rqc"))
install.packages(c("parallel"))

Introduction

Comment le séquencage est effectué?

Crédit: J. Couillard

DOI: 10.1016/B978-0-12-802234-4.00002-1

Comment le séquencage est effectué?

Note

  • Chaque nucléotide à un fluorophore spécifique (4 canaux).
  • Certains appareils (ex. NextSeq) vont utiliser deux canaux car c’est plus rapide et moins coûteux.

DOI: 10.1016/B978-0-12-802234-4.00002-1

Le signal renvoyé par le laser est traduit en chromatographe pour chaque amplicon séquencé.

Comment le séquencage est effectué?

Exemple de chromatographe

Source: LabXchange.org

Comment le séquencage est effectué?

Exemple de chromatographe

Important

Le score est basé sur:

  • La séparation entre les clusters (pureté du signal)
  • Le rapport signal/bruit
  • La distribution de l’intensités des 4 canaux
  • La position dans le read (la qualité diminue généralement vers la fin)

Source: wikipedia.org

Fichiers de sortie du séquenceur Illumina

Exemple de nomenclature de fichier de sortie du séquenceur Illumina: sample1_1_L001_R1_001.fastq.gz

Signification Exemple
sample1 Identifiant de l’échantillon sample1, sample2
_1 Numéro de réplicat 1, 2, 3…
L001 Numéro de ligne / lane: L001-L008
R1/R2 Direction de lecture R1=avant, R2=arrière
001 Segment de fichier 001, 002…

Note

  • Séquençage paired-end = fichiers R1 + R2 pour chaque échantillon avec l’ensemble des amplicons.
  • Numéro de Ligne ou lane = Une flowcell Illumina est divisée en plusieurs voies physiques parallèles où le séquençage se déroule simultanément.

Structure du fichier FASTQ

Contient plusieurs lignes dont un ensemble de 4 lignes correspond à un amplicon.

system.file(package="ShortRead", "extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz") |> 
  gzfile() |> 
  readLines(n = 12)
 [1] "@ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1"                  
 [2] "GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGACATGAAGTCAATGAAGGCCTGGAATGTCACTACCCCCAG"
 [3] "+"                                                                       
 [4] "HHHHHHHHHHHHHHHHHHHHEBDBB?B:BBGG<DDAA?AABFEFBDBD@DDECEE3>:?;@@@>?=BAB?##"
 [5] "@ERR127302.21406531 HWI-EAS350_0441:1:88:9330:2587#0/1"                  
 [6] "CTAGGGCAATCTTTGCAGCAATGAATGCCAATGGGTAGCCAGTGGCTTTTGAGGCCAGAGCAGACCTTCGGG"
 [7] "+"                                                                       
 [8] "IIIIHIIIGIIIIIIIHIIIIEGBGHIIIIHGIIHIIIIIIIHIIIHIIIIIGIIIEGIIGBGE@DDGGGIG"
 [9] "@ERR127302.22173106 HWI-EAS350_0441:1:91:10434:14757#0/1"                
[10] "TGGGCTGTTCCTTCTCACTGTGGCCTGACTAAAACCCAGTGGCATTAAGAAAGAGTCACGTTTCCCAAGTCT"
[11] "+"                                                                       
[12] "GGHBHGBGGGHHHHDHHHHHHHHHFGHHHHHHHHHHHHHHHHHHHHHGHFHHHHHHHHHHHHHH8AGDGGG>"

Structure du fichier FASTQ

system.file(package="ShortRead", "extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz") |> 
  gzfile() |> 
  readLines(n = 4)
[1] "@ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1"                  
[2] "GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGACATGAAGTCAATGAAGGCCTGGAATGTCACTACCCCCAG"
[3] "+"                                                                       
[4] "HHHHHHHHHHHHHHHHHHHHEBDBB?B:BBGG<DDAA?AABFEFBDBD@DDECEE3>:?;@@@>?=BAB?##"
  1. Ligne [1] : Identifiant de séquence (commence par @)
  2. Ligne [2] : Séquence d’ADN (A, T, G, C) de l’amplicon
  3. Ligne [3] : Séparateur (+)
  4. Ligne [4] : Scores de qualité (un par base)

Ligne d’identification

Décomposition de l’identifiant

@SH00321:6:BWR98207-2813:1:1101:1065:1015 1:N:0:AACCATAGAA+GGCGAGATGG
  • SH00321 : ID de l’instrument
  • 6 : Numéro d’exécution
  • BWR98207-2813 : ID de la flowcell
  • 1:1101:1065:1015 : Tuile et coordonnées X:Y
  • 1:N:0 : Numéro de lecture, flag de filtre, bits de contrôle
  • AACCATAGAA+GGCGAGATGG : Codes-barres d’échantillon (double indexation)

Interprétation des scores de qualité

Exemple de chaine de caractères donnant le score de qualité (1 caractère par nucléotide): GGGGGGGGGGGGGGGGGGGGGGGG9GG-G9G9G9GGGG

Caractère ASCII Score Phred Taux d’erreur Précision
G 38 0,016% 99,984%
9 24 0,40% 99,60%
- 12 6,31% 93,69%

Tip

Règle générale : Q30 ou plus = haute qualité (99,9% de précision)

Scores de qualité du séquencage

Les scores de qualité utilisent l’encodage ASCII.

 Quality encoding: !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHI
                   |         |         |         |         |
    Quality score: 0........10........20........30........40

Scores de qualité du séquencage

Avec R, comment convertir l’encodage ASCII en score de Phred et en probabilité d’obtenir une erreur (précision)?

# Convertir la chaîne de qualité en scores Phred
quality_string <- "GGGGGGG99GG-GGG"
quality_scores <- as.integer(charToRaw(quality_string)) - 33

# Afficher la correspondance
data.frame(
  Caractere = strsplit(quality_string, "")[[1]],
  Score_Phred = quality_scores,
  Precision = paste0(round((10^(-quality_scores/10)) * 100, 2), "%")
)

Scores de qualité du séquencage

   Caractere Score_Phred Precision
1          G          38     0.02%
2          G          38     0.02%
3          G          38     0.02%
4          G          38     0.02%
5          G          38     0.02%
6          G          38     0.02%
7          G          38     0.02%
8          9          24      0.4%
9          9          24      0.4%
10         G          38     0.02%
11         G          38     0.02%
12         -          12     6.31%
13         G          38     0.02%
14         G          38     0.02%
15         G          38     0.02%

Lire un fichier FASTQ avec R

# Lire le fichier FASTQ
fq <- system.file(package="ShortRead", "extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz") |> ShortRead::readFastq()

# Examiner la structure
fq
class: ShortReadQ
length: 20000 reads; width: 72 cycles

Extraire les informations principales

# Nombre d'amplicons
length(fq)
[1] 20000
# Extraire les amplicons
sequences <- ShortRead::sread(fq)

# Extraire les identifiants
ids <- ShortRead::id(fq)

# Extraire les scores de qualité
qualities <- Biostrings::quality(fq)

# Premier identifiant et séquence
ids[1]
BStringSet object of length 1:
    width seq
[1]    53 ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1
sequences[1]
DNAStringSet object of length 1:
    width seq
[1]    72 GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGAC...GTCAATGAAGGCCTGGAATGTCACTACCCCCAG

Résumé statistique d’un fichier FASTQ

# Longueur des séquences
seq_lengths <- Biostrings::width(sequences)
summary(seq_lengths)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     72      72      72      72      72      72 
# Distribution des scores de qualité
qual_matrix <- as(qualities, "matrix")
mean_quality <- rowMeans(qual_matrix)
summary(mean_quality)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.00   33.53   37.35   34.84   38.85   39.99 
# Composition en bases
Biostrings::alphabetFrequency(sequences, as.prob = TRUE)[1:4, 1:4]
             A         C         G         T
[1,] 0.1944444 0.2638889 0.2916667 0.2500000
[2,] 0.2361111 0.2222222 0.3194444 0.2222222
[3,] 0.2361111 0.2638889 0.2222222 0.2777778
[4,] 0.2500000 0.2638889 0.1666667 0.3194444

Assurance qualité avec la librairie Rqc

fq_files <- system.file(package="ShortRead", "extdata/E-MTAB-1147/") |> list.files(full.names = TRUE)
qa <- Rqc::rqcQA(fq_files)
Rqc::rqcCycleQualityPlot(qa)

Analyser la qualité du séquencage

Pourquoi la qualité des lectures décroit au fur et à mesure des itérations?

  1. Dégradation des fluorophores
  2. Accumulation d’erreurs de phase, soit la désynchronisation des clusters
  3. Épuisement des réactifs, compétition pour les nucléotides par exemple
  4. Dommages à l’ADN par stress chimique par exemple
  5. Accumulation de sous-produits (ex. Résidus de réaction s’accumulent)

Pratique - Analysez vos propres données

Note

Exercice simple: Utilisez le code R que nous venons de voir pour analyser vos propres fichiers FASTQ

  1. Créer un dossier d’analyse
  2. Placez vos fichiers FASTQ dans sous-dossier data/
  3. Exécutez les chunks (le code des diapos) pour analyser la qualité
  4. Comparez les résultats avec les exemples montrés
  5. Bonus: Explorer les autres graphiques dans le package Rqc

Tip

Si vous n’avez pas de données, vous pouvez télécharger l’exemple ici. Ce sont 10 échantillons de métabarcoding mitofish-12S

10:00

Points clés à retenir

  • FASTQ = format standard pour les données de séquençage
  • 4 lignes par lecture d’amplicon : ID, lecture de l’amplicon, séparateur, qualité de la lecture
  • Paired-end : les fichiers R1 et R2 contiennent les lecture d’amplicons à apparier
  • Scores de qualité : Plus élevé = meilleur (viser Q30+)
  • La qualité se dégrade vers la fin des lectures (normal)

Qu’est ce que BARQUE?

Un large éventail d’outils en bioinformatique

Hakimzadeh et al, 2024

Qu’est ce que BARQUE?

  • Programme écrit par Éric Normandeau (IBIS, ULaval)
  • Disponible sur Github depuis 2017
  • C’est un pipeline de traitement de données de séquences spécialisé sur le metabarcoding
  • Il agit comme chef d’orchestre d’un ensemble de programmes

Qu’est ce que BARQUE?

Creative Commons
  • FLASH v1.2.11+: fusion des lectures paired-end
  • VSEARCH v2.14.2+: outil d’analyse de séquences (version obligatoire)
  • Trimmomatic: outil de filtrage et découpage des lectures
  • NCBI BLAST: annotation des séquences non-correspondantes
  • Un ensemble de programmes UNIX comme cut, sort, uniq etc. qui manipulent des chaînes de caractères ou des nombres
  • R & Python: Visualisation et synthèse des résultats

Les grandes étapes de BARQUE

Vue d’ensemble du pipeline de métabarcoding

Étape 1: Validation du projet

  • Vérification des fichiers d’entrée
  • Validation de la nomenclature : SampleID_*_R1_001.fastq.gz
  • Vérification des logiciels
  • Présence de la base de données correspondante à l’amorce

Étape 2: Nettoyage des données

  • Nettoyage des lectures avec Trimmomatic
  • Suppression des bases de faible qualité
  • Découpage selon les paramètres définis
  • Retrait des amplicons non-appariés

Étape 3: Appariement des lectures

  • Appariement des amplicons (R1 + R2) avec FLASH
  • Création d’une lecture consensus pour chaque amplicon
  • Amélioration de la qualité globale

Étape 4: Séparation des amplicons

Pour chaque lecture fusionnée:

  • Recherche du primer forward au début
  • Recherche du primer reverse à la fin (complément inverse)
  • Calcul de la distance avec les primers de référence
  • Vérification de la taille d’amplicon attendue

Étape 5: Retrait des chimères

  • Détection des séquences chimériques avec VSEARCH
  • Les chimères sont des artefacts de PCR
  • Conservation des séquences authentiques

Étape 6: Assignement taxonomique

  • Attribution des espèces via bases de données correspondantes à l’amorce
  • Utilisation de BOLD (COI), MitoFish (12S), etc.
  • Génération des tables d’occurrence par espèce

Télécharger BARQUE

Pratique

Télécharger BARQUE en cliquant sur le lien

Ou en rendant à l’adresse suivante: https://github.com/enormandeau/barque

Les étapes en tant qu’utilisateur

  1. On place les fichiers fastq.gz avec la bonne nomenclature dans le dossier 04_data
  2. On configure les amorces (P3 et P5) utilisées avec le fichier 02_info/primers.csv
  3. On ajuste les paramètres du pipeline avec le fichier 02_info/barque_config.sh
  4. On s’assure que la base de données BOLD est téléchargé si on travaille avec les amorces COI.
  5. On éxecute le pipeline

Tip

Si vous voulez réutiliser la configuration du pipeline pour d’autres fois, c’est une bonne pratique de dupliquer le fichier et de le sauvegarder sous un autre nom.

Étape 1 - Placer les fichiers dans 04_data

On peut faire tout de suite l’étape 1 avec vos données.

Tip

Si vous avez pas de données, vous pouvez télécharger celles fournies par Étienne Normandeau, en cliquant ici. Ce sont 10 échantillons de métabarcoding mitofish-12S, chacun.

Étape 2 - Configuration des amorces

On configure les amorces (P3 et P5)

  • Ouvrir le fichier 02_info/primers.csv
  • Retirer le symbole # pour la ligne correspondante à votre amorce
  • Laisser le symbole # sur les autres lignes

Étape 3 - Ajustements des paramètres du pipeline

# Modify the following parameter values according to your experiment
# Do not modify the parameter names or remove parameters
# Do not add spaces around the equal (=) sign
# It is a good idea to try to run Barque with different parameters 

# Global parameters
NCPUS=20                    # Number of CPUs to use. A lot of the steps are parallelized (int, 1+)
PRIMER_FILE="02_info/primers.csv" # File with PCR primers information

# Skip data preparation and rerun only from vsearchp
SKIP_DATA_PREP=0            # 1 to skip data preparation steps, 0 to run full pipeline (recommended)

# Filtering with Trimmomatic
CROP_LENGTH=200             # Cut reads to this length after filtering. Just under amplicon length

# Merging reads with flash
MIN_OVERLAP=30              # Minimum number of overlapping nucleotides to merge reads (int, 1+)
MAX_OVERLAP=280             # Maximum number of overlapping nucleotides to merge reads (int, 1+)

# Extracting barcodes
MAX_PRIMER_DIFF=8           # Maximum number of differences allowed between primer and sequence (int, 0+)

# Running or skipping chimera detection
SKIP_CHIMERA_DETECTION=0    # 0 to search for chimeras (RECOMMENDED), 1 to skip chimera detection
                            #   or use already created chimera cleaned files

# vsearch
MAX_ACCEPTS=20              # Accept at most this number of sequences before stoping search (int, 1+)
MAX_REJECTS=20              # Reject at most this number of sequences before stoping search (int, 1+)
QUERY_COV=0.6               # At least that proportion of the sequence must match the database (float, 0-1)
MIN_HIT_LENGTH=100          # Minimum vsearch hit length to keep in results (int, 1+)

# Filters
MIN_HITS_SAMPLE=10          # Minimum number of hits in at least one sample to keep in results (int, 1+)
MIN_HITS_EXPERIMENT=20      # Minimum number of hits in whole experiment to keep in results (int, 1+)

# Non-annotated reads
NUM_NON_ANNOTATED_SEQ=200   # Number of unique most-frequent non-annotated reads to keep (int, 1+)

# Multiple hits
MIN_DEPTH_MULTI=10          # Min depth to report unique reads per sample in multiple hit reports

# OTUs
SKIP_OTUS=1                 # 1 to skip OTU creation, 0 to use it
MIN_SIZE_FOR_OTU=20         # Only unique reads with at least this coverage will be used for OTUs

On ajuste les paramètres du pipeline dans le fichier 02_info/barque_config.sh

Tip

Pour déterminer le nombre de coeur sur votre ordinateur, utiliser la commande R suivante parallel::detectCores() - 2

Étape 4 - Base de données par type d’amorces

BARQUE utilise différentes bases de données selon le type d’amorces utilisé pour l’amplification:

Bases de données intégrées:

  • MitoFish (12S rRNA): Séquences mitochondriales de poissons
  • BOLD (COI): Base de données génétiques du code-barres de la vie
  • SILVA (18S rRNA): Séquences ribosomales eucaryotes
  • Bases personnalisées: Possibilité d’ajouter vos propres bases

Formats requis:

  • Fichiers FASTA avec headers normalisés
  • Taxonomie incluse dans le header
  • Placement dans 03_databases/

Étape 4 - Télécharger des base de données

BOLD (COI) et SILVA (18S)

Les bases de données BOLD et SILVA pré-formatées pour BARQUE sont disponibles ici:

https://www.ibis.ulaval.ca/services/bioinformatique/barque_databases/

Important

Après téléchargement, renommer le fichier en silva.fasta.gz ou ****.fasta.gz** et le placer dans 03_databases/

On est prêt à executer BARQUE !!

Pourquoi utiliser docker

BARQUE requière plusieurs programmes qui ne sont pas tous compatibles avec l’environnement Windows

Hakimzadeh et al, 2024

Comment docker fonctionne?

Docker crée un environnement isolé (conteneur) qui contient:

  • Tous les programmes requis par BARQUE
  • Les bonnes versions des dépendances
  • Un système Linux compatible

Le conteneur peut accéder aux fichiers de votre ordinateur via des volumes montés.

Tip

Docker permet d’utiliser BARQUE sur Windows, Mac ou Linux sans installer manuellement tous les programmes requis

Comment docker fonctionne?

@Soroosh Nazem

La recette

Déclaration de l’image docker

# Base image with Ubuntu
FROM ubuntu:24.04

# Update and install essential packages
RUN apt-get update && apt-get install -y \
    wget \
    curl \
    build-essential \
    software-properties-common \
    default-jre \
    parallel \
    r-base \
    bc \
    python3-setuptools \
    ca-certificates \
    && rm -rf /var/lib/apt/lists/*

# Install FLASH (v1.2.11+)
RUN wget https://github.com/dstreett/FLASH2/archive/refs/tags/2.2.00.tar.gz && \
    tar -xzf 2.2.00.tar.gz && \
    cd FLASH2-2.2.00 && \
    make && \
    mv flash2 flash && \
    cp flash /usr/local/bin && \
    cd .. && rm -rf FLASH2-2.2.00 2.2.00.tar.gz

# Install VSEARCH v2.14.2+
RUN wget https://github.com/torognes/vsearch/releases/download/v2.30.0/vsearch-2.30.0-linux-x86_64.tar.gz && \
    tar -xzf vsearch-2.30.0-linux-x86_64.tar.gz && \
    cp vsearch-2.30.0-linux-x86_64/bin/vsearch /usr/local/bin && \
    rm -rf vsearch-2.30.0-linux-x86_64*

# Default shell
CMD ["/bin/bash"]

Étape 1: Télécharger l’image Docker

Ouvrir un terminal (MacOS) ou PowerShell (Windows) exécuter la commande suivante pour télécharger l’image BARQUE:

docker pull steveviss/barque:latest

Cette commande télécharge l’image Docker contenant BARQUE et tous ses programmes requis (FLASH, VSEARCH, Trimmomatic, etc.)

Note

Le téléchargement peut prendre quelques minutes selon votre connexion Internet

Étape 2: Exécuter l’image Docker en mode interactif

Pour lancer BARQUE dans un conteneur Docker avec accès à vos fichiers:

# Pour macOS
docker run -it -v /chemin/vers/barque:/barque steveviss/barque:latest

# Pour Windows
docker run -it -v C:\chemin\vers\barque:/barque steveviss/barque:latest

Important

Remplacer /chemin/vers/barque par le chemin absolu vers votre dossier BARQUE local

Explication des options:

  • -it : Mode interactif avec terminal
  • -v /chemin/vers/barque:/barque : Monte votre dossier local dans le conteneur
  • steveviss/barque:latest : L’image Docker à utiliser

Étape 3: Exécuter BARQUE dans le conteneur Docker

Une fois dans le conteneur Docker, vous êtes dans un environnement Linux avec BARQUE prêt à l’emploi:

# Naviguer vers le dossier BARQUE
cd /barque

# Lancer BARQUE avec votre fichier de configuration
./barque 02_info/barque_config.sh

Pour quitter le conteneur:

exit

Application Shiny

Exécuter BARQUE de manière embarqué

Application disponible à l’adresse suivante: https://cloud.taq.info, contacter moi pour un accès: steve.vissault@inrs.ca

Application Shiny

Exécuter BARQUE de manière embarqué

Application disponible à l’adresse suivante: https://cloud.taq.info, contacter moi pour un accès: steve.vissault@inrs.ca

Rapport de metabarcoding

Obtenir le rapport sur les résultats de BARQUE

Rapport de metabarcoding

Installation de barqueReport

Vous pouvez installer la version de développement de barqueReport depuis GitHub avec :

install.packages("remotes")
remotes::install_github("taq-community/barqueReport")

Rapport de metabarcoding

Générer un rapport de métabarcoding

La fonction principale génère un rapport HTML interactif à partir des résultats du pipeline BARQUE:

library(barqueReport)

# Sélection interactive du dossier BARQUE et du fichier de qualité Illumina
generate_metabarcoding_report(
  barque_output_folder = file.choose(),  # Sélectionner le dossier de sortie barque
  samples_ids = c("ST1", "ST2", "ST3", "ST4"),
  title = "Étude du Lac Supérieur",
  subtitle = "",
  author = "Steve Vissault, Julie Couillard & Tuan Ahn To",
  illumina_quality_file = file.choose(),  # Sélectionner le fichier Quality_Metrics.csv
  blank_lab = "LAB_BLANK",
  blank_field = "FIELD_BLANK"
)